Morphogen gradient from a noisy source 
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We investigate the effect of time-dependent noise on the shape of a morphogen gradient in a 
developing embryo. Perturbation theory is used to calculate the deviations from deterministic 
behavior in a simple reaction-diffusion model of robust gradient formation, and the results are 
confirmed by numerical simulation. It is shown that such deviations can disrupt robustness for 
sufficiently high noise levels, and the implications of these findings for more complex models of 
gradient-shaping pathways are discussed. 



Successful embryonic development requires that dif- 
ferent regions within the embryo receive distinct sets of 
developmental instructions. These instructions are of- 
ten relayed by a morphogen, a signaling molecule that 
induces different cell fates at different concentrations. 
By establishing concentration gradients of various mor- 
phogens during early stages of development, the embryo 
provides each cell with the positional information it needs 
for the proper implementation of a body plan (Fig. 

III 

A morphogen gradient must reliably demarcate pre- 
cise and accurate boundaries between groups of cells de- 
spite substantial genetic and environmental variability. 
Past studies have argued that this task might be accom- 
plished in at least two different ways. First, it has been 
proposed that the downstream pathway that responds 
to a morphogen may have built-in mechanisms for in- 
terpreting sharply-defined positional information from a 
gradient that varies quite substantially from embryo to 
embryo 0. Second, the shape of a gradient might it- 
self be "robust" to certain kinds of interference. Eldar et 
al. have presented evidence that the shape of the BMP 
morphogen gradient may not be altered significantly by 
several types of mutations in the pathway responsible 
for establishing the gradient 3]. Furthermore, their the- 
oretical investigations of dynamical models for several 
gradient-shaping pathways have uncovered steady-state 
solutions that are effectively independent of certain ki- 
netic parameters and initial or boundary conditions 0,01 • 
Thus, there is some indication that a pathway's architec- 
ture could ensure that the shape of the resulting gradient 
did not depend on the rates of synthesis or the concen- 
trations of some or all components of the pathway. 

The theoretical studies of robustness mentioned above 
rely on deterministic reaction-diffusion equations that 
treat the influxes of the components of a gradient-shaping 
pathway as constant, time-independent quantities. How- 
ever, a growing number of experimental studies have 
demonstrated that gene expression at the cellular level 
is an inherently "noisy" process, during which the rates 
of synthesis of various gene products fluctuate substan- 



FIG. 1: Morphogen concentration gradients cj>o(x) are plotted 
from © for two different values of the production rate r/o. 
Where <j)o is greater than the threshold concentration <j>* , a 
cell fate of type I is induced (unshaded). In contrast, <j>o < 
4>* induces a cell fate of type II (shaded). The process by 
which the gradient is shaped is said to be robust to rjo because 
both curves demarcate essentially the same location x* for the 
boundary between cell types. 



tially about their average values over time 0, 0, S @ • I n 
the presence of such noise, a gradient would be constantly 
driven away from its steady-state shape, and might even 
adopt a different shape, on average. Thus, the question 
of whether deterministic models can adequately capture 
the dynamics of morphogens that are, in reality, pro- 
duced and secreted stochastically is an open one whose 
answer may bear profoundly on our understanding of how 
a viable final product is assured during the process of em- 
bryonic development. 

This Letter introduces the novel dimension of time- 
dependent noise to the study of morphogen gradients. 
A perturbative treatment of a simple reaction-diffusion 
model for gradient formation is used to calculate both 
the noise-averaged shift of the system away from its "ro- 
bust" deterministic steady-state, and the system's fluc- 
tuations about this average. These analytical results are 



subsequently confirmed by a numerical simulation. Com- 
paring our findings to the steady-state solution of the 
reaction-diffusion model, we find that deterministic ro- 
bustness can indeed be compromised by noise, but only 
when the intensity of noise is sufficiently high. Finally, we 
argue that models of robust gradient formation that in- 
corporate a more realistic number of components should 
be even more susceptible to the effects of noise than the 
simple model investigated here. 

We begin by considering a single morphogen of local 
concentration <f>(x, t) that is produced at the origin at 
a constant rate r) and diffuses out to infinity along the 
positive x-axis with Fick's constant D. The morphogen 
accelerates its own rate of degradation such that it is re- 
moved from the system at a local rate of f4>{x, t) 2 , which 
gives the equation of motion 



d t <P = D dl<j> - fcj? + m 5(x) 



(1) 



Experimental evidence suggests that this mechanism of 
self-enhanced ligand degradation is a qualitatively cor- 
rect description of how the Wg gradient is established in 
fruit flies Q. 

When the gradient reaches a steady-state shape (j>o, it 
will divide space into two regions: one where the concen- 
tration of morphogen is higher than some threshold value 
(f>* and a cell fate of type 1 is induced, and one where the 
concentration is lower and a cell fate of type 11 is induced 
(Fig. 2]). The boundary between the two regions lies at 
the coordinate x* that satisfies (fio(x*) = <j>* . 

The steady-state solution that satisfies d t 4>o = is 
given by 
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The boundary coordinate is 
'(0)- 



therefore given by x* = y/6D/f(j)* — e = x* (0) - e. If e < 
x*(e), then <ilogx*/dloge will be very small, meaning 
that relatively large percentage changes in e will only 
introduce minor percentage shifts in x* . Thus, so long as 
e -C x*(0), the location of a cell type boundary specified 
as in Figure Q] will be robust to the concerted shifts in i]q 
that might result from genetic or environmental variation 
within a population of embryos Q . 

The preceding analysis operated under the assump- 
tion that the rate at which the morphogen is synthesized 
and secreted is constant over time in a single embryo. 
We now relax this assumption in order to study the con- 
sequences of introducing time-dependent noise into the 
system. The first step in doing this is to define new quan- 
tities <j>(x,t) = 4>o{x) +S(f>(x,t) and r)(t) = rj +Si](t). A 
rescaling of space, time, noise, and concentration in the 
interest of consolidating free parameters results in the 
new equation of motion 

d t S(j) = d 2 Jct) - - 12 .. <ty - 8(j) 2 + Srj S(x) (3) 
(x + ey 




FIG. 2: The Feynman diagram for a general, nth-order pro- 
cess, with time increasing from left to right, n distinct fluctua- 
tion events at the origin propagate through space and interact 
with each other to produce a positive contribution to A at x' 
(inner box). The self-interaction of this variance summed over 
all points x' then induces a non-zero {S(f>) at x (outer box). 



We have computed the mean shift (5(f>(x)) and vari- 
ance A = (8<fi(x) 2 ) — (S(j)(x)) 2 (Here, (. . .) denotes aver- 
aging over all realizations of the noise Sr](t)). Since di- 
mensional analysis requires that these two quantities de- 
pend on e(ryo), they are a means to gauge the impact that 
noise has on the robustness of the gradient to changes in 
r)Q. The noise itself is taken to satisfy (Sr)(t)) — and 
(5rj(t)Srj(t')} — 2^S(t — t'). A power-counting argument 
shows that the effects of more realistic noise correlations 
should not alter the behavior for large x. 

The quantities of interest can be calculated using per- 
turbation theory. The first step is to construct a Green's 
propagator, which is composed of the eigenfunctions 'ipx 
of the linear operator — d 2 + 12/(x + e) 2 in the right-hand 
side of 



G(xi,x 2 ;uj) 



i^\(xi)^\(x 2 ) 
—iuj + A 



(4) 



The eigenfunctions may be approximated as ip\(x) = 
\J (x + e)/2 J7/2 [V\(x + e)]. While terms proportional 
to Bessel functions of the second kind are needed for ex- 
actness when e ^ 0, their contribution should be small 
for the soft modes of the system that contribute most 
strongly to averages over the noise, particularly in the 
case of interest where e is small compared with the size 
of the system. 

The iterative solution to J3J obtained by repeated sub- 
stitution for the non-linear term can be expressed as a 
sum over tree Feynman diagrams rooted at x, each edge 
corresponding to G, each node to the non-linear inter- 
action, and the end of each branch corresponding to an 
insertion of the noise. Averaging over the noise then has 



the effect of pairing these ends in all possible ways. The 
nth-order contribution to A is represented in the inner 
box of Figure [21 The two lines that converge on x' must 
originate at vertices located at points x\ and x%, which 



can be assumed to be less than x' for reasons that will 
be explained below. Integration over to and the propa- 
gator eigenvalue for G(xi,x') reveals that the nth-order 
diagram contains the integral 



J d\ {x 1 + e) y /{x 1 +e){x 2 +e)J 7/2 {VX[x 1 + e])J 7/2 {V\[x' + e])/ 7/2 (\/A> 2 + e])K 7/2 {V\[x' + e]) (5) 

I 



Assuming that x' +e 3> Xi+e for all points Xi that tempo- 
rally precede x' in Fig. |2 the heavy damping from K 7 i 2 
permits the substitution J 7 / 2 (\f~\[xi + e])I 7 / 2 {\f~\[x 2 + 
e]) ~ X 1 / 2 . The subsequent integral over A in comes 
out the same regardless of the order in perturbation the- 
ory. Thus, we obtain A ~ l/(x') 8 . 

The above result was based on the assumption that x' 
is always much larger than the spatial coordinates of ver- 
tices located at earlier time coordinates. Obviously, when 
one integrates over all possible vertex positions, one will 
be forced to consider cases where x\ 3> x' . However, 
free fluctuations constantly relax towards 5<j) = as they 
propagate through space, and their effects are clearly the 
strongest if they interact non-linearly prior to undergo- 
ing most of this relaxation. Thus, by assuming that all 
interactions contributing to a diagram occur before prop- 
agation over any significant distance takes place, we ig- 
nore processes that would have contributed only weakly 
to the final answer. 

With A in hand, we now proceed to the outer box of 
Fig. [3 and calculate the average shift in 5(f> using 

{5<f>(x)) = f dx'A(x')G(x,x';u = 0) (6) 

Since 



G(x>x';0) = ^(^j ~ (7) 

the leading term in becomes (5<f>(x)) ~ 1/x 3 . It 
should be noted, however, that the cubic power law is 
not a result of the specific form of the fluctuations A as 
long as they fall off sufficiently rapidly with x' . Rather, 
the observed decay would be brought in as part of any 
computation of JSJl because of the cubic spatial depen- 
dence of G(x > x';0). Thus, it seems likely that little 
was sacrificed in generality by assuming a white spec- 
trum for the noise. 

The system defined in llll.'il) was simulated in order 
to test whether {5<j>) and A would exhibit the predicted 
power law behaviors. As Figure [3] illustrates, the con- 
currence of the data with the analytical calculations is 
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FIG. 3: The equation of motion was solved numerically 
for At = .001, Ax = .2, / = 5, 7 = .2, and D — 1 on 
< x < 8.8. Boundary conditions were imposed such that 
9x5(j>\ x _ = 6(f>(8.8) = 0. Here, we plot the inverse cube root 
of the magnitude of the average shift in 5(j>, along with the 
inverse eighth root of the fluctuation of the shift about its av- 
erage value (inset). Regression analysis of the solid lines gave 
linear correlation coefficients > .998; the dotted lines repre- 
sent regions where the finite-size effect due to the boundary 
condition at x — 8.8 becomes noticeable. 

essentially exact, except near x = 8.8, where a bound- 
ary condition requiring 8<j) to vanish artificially skews the 
curve towards zero. 

We are now in the position to assess whether the 
gradient remains robust to changes in the average pro- 
duction rate (77) = ??o in the presence of white noise of 
strength 7. Recall that the condition for robustness to 770 
in the steady-state case was e <§; x*(0). Now examining 
the noise-averaged shape of the gradient, we find that 

. , . ., 6D/f A 
0(a) ~ - — (8) 

(x + e) z x a 

where A(rjo, 7) is an unknown, positive amplitude that in- 
creases in a monotonic, unbounded fashion with 7 from 
the value of zero at 7 = 0. Comparing to J5J, it is ap- 
parent the impact of the noise at leading order is that 
it drives up the effective value of e. This means that 
although our criterion for robustness may be met in the 
absence of noise, there must be a 7* (770) for which 7 > 7* 



implies that 2e(r] ) + A(j, r/ Q ) > x*(0). Thus, the robust- 
ness of the system to changes in t/q is reduced in the 
presence of noise, and can be eliminated altogether if the 
noise level is sufficiently high. 

It is also interesting to consider the effects that the 
fluctuations might themselves have on the process of de- 
velopment. If the embryo had an infinite amount of time 
to sample the fluctuating level of morphogen at each 
point in space, it could "measure" and respond to the 
gradient {4>) with perfect precision. Instead, positional 
information must be recovered from a morphogen gradi- 
ent over the course of some finite period T dictated by 
the timetable of embryonic growth. A calculation simi- 
lar to the one carried out in J5J reveals that the variance 
measured at x during this period will be less than A and 
asymptotically proportional to T~ 1 x~ 6 (The slower x~ 6 
decay of fluctuations results from the fact that points 
more distant from the origin receive information about 
the effective value of e more slowly). Thus, the precision 
of cell type boundaries must be limited approximately by 
Aa;* ~ (x*)~ 3 /((f>'(x*)). This uncertainty in the bound- 
ary position should grow monotonically with 7 and de- 
pend only very weakly on x* . Thus, regardless of where 
the boundary lies on average, a sufficiently high noise 
strength could push the precision in x* below tolerable 
levels. 

The above discussion focused on a reduced model 
of morphogen gradient formation that involves only one 
chemical species. In fact, even the simplest realizations 
of self-enhanced ligand degradation in vivo involve in- 
termediary receptors and proteases that enable the mor- 
phogen to regulate its own rate of degradation 4] . How- 
ever, there are at least two reasons to believe that this 
additional sophistication should only magnify the effects 
we have observed in the reduced model. Firstly, this 
study has focused on a system in which the only source 
of noise was localized at the single point in space at which 
the damping of fluctuations was strongest. In contrast, 
more realistic models may contain multiple sources of 
noise that extend across the entire system. Secondly, in 
more complex pathways, it is possible for the concen- 
tration of one chemical species to control the relaxation 
time of fluctuations in a second species. Thus, a path- 
way that is robust to the concentration of a component 
in the deterministic steady-state may exhibit fluctuations 
and an induced shift away from the steady-state that are 
not robust to this concentration. For example, we have 
recently shown that under reasonable assumptions, the 
gradient shape established by the BMP pathway stud- 
ied in |3] should exhibit increased variability when the 
ambient concentration of the short gastrulation protein 
is reduced. This predicted variability is anecdotally con- 



firmed by experiments carried out in [9|, but the issue 
still demands a more quantitative investigation. 

This study is, to our knowledge, the first to eluci- 
date the impact that time-varying stochasticity in gene 
expression can have on the shaping of morphogen gradi- 
ents. Using both analytical and numerical techniques, we 
have shown that noise can interfere with the putatively 
robust specification of positional information by induc- 
ing a non-robust mean shift in the gradient away from 
its steady-state shape, and can also cause fluctuations in 
the position of a cell fate boundary. There is furthermore 
good reason to believe that these effects should be more 
pronounced in models more realistic than the simplified 
system considered here. However, while it has been made 
clear in this work that noise has potential significance, 
it remains a question for experimenters whether fluctu- 
ations are large enough that they can play a detectable 
role in the formation of morphogen gradients. Yet even in 
the cases where noise turns out to have little effect, there 
will still remain the intriguing possibility that specific 
mechanisms are responsible for attenuating fluctuations 
in order to prevent them from interfering substantially 
with the action of morphogens. In either event, our find- 
ings should provide motivation for an exciting new line of 
experimental inquiry that may contribute to a richer un- 
derstanding of the factors that play a decisive role during 
embryonic development. 
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